#ifndef FJF_NN_SSE_COR_N_01_H
#define FJF_NN_SSE_COR_N_01_H

#include <iostream>
#include <math.h>
#include <xmmintrin.h>
#include <stdio.h>
#include <string.h>

//my includes
#include "utils.02.h"

namespace dj01 {
	
class FJF_NN2_SSE_COR_n {	
	
	public:	
	
	// n, j1 should be even, j2 should be odd, arra
	inline virtual void ff(	uint n, uint N, 
									uint i1_,uint j1_,uint i2_,uint j2_, 
									double*m, 
									double*x0,double*x1,double*x2,
									double*v0,double*v1,double*v2,  
									double*f0,double*f1,double*f2, 
									double*j0,double*j1,double*j2, 
									double*pot, double eps2) {

		if (n%2!=0 && n>1) n++; // we need even number for SSE
		
		uint dj=10000; // array indexing range (must be align 2) 10000 is optimal
		uint kn = ceil((double)N/dj);
		uint _i1, _i2, _j1, _j2, _j11, _j21, _j12, _j22, jL, jR ;	
	
		double tmp[2] __attribute__((aligned(16)));	

		__m128d _eps2 =_mm_setr_pd(eps2, eps2);
		__m128d _3 =_mm_setr_pd(3, 3);
			
		for (uint i=0; i<=N; i++)
			f0[i]=f1[i]=f2[i]=j0[i]=j1[i]=j2[i]=pot[i]=0;
		
		for (uint k=0; k<kn; k++) {  // lets devide j walk on kn parts with size dj performance reasons

			jL= k*dj; // lets define Light j-coordinate 
			N-jL>dj ? jR=jL+dj-1 : jR = N-1; // j-Left coordinate
				
			_i1=i1_; _i2=i2_; // ??? should be modified later on  ???
			
			for (uint i=_i1; i<=_i2; i++) {
				
				__m128d mi =_mm_setr_pd(m[i], m[i]);
				
				__m128d x0i=_mm_setr_pd(x0[i], x0[i]);
				__m128d x1i=_mm_setr_pd(x1[i], x1[i]);
				__m128d x2i=_mm_setr_pd(x2[i], x2[i]);
			
				__m128d v0i=_mm_setr_pd(v0[i], v0[i]);
				__m128d v1i=_mm_setr_pd(v1[i], v1[i]);
				__m128d v2i=_mm_setr_pd(v2[i], v2[i]);
			
				__m128d f0i=_mm_setr_pd(0, 0);
				__m128d f1i=_mm_setr_pd(0, 0);
				__m128d f2i=_mm_setr_pd(0, 0);
			
				__m128d j0i=_mm_setr_pd(0, 0);
				__m128d j1i=_mm_setr_pd(0, 0);
				__m128d j2i=_mm_setr_pd(0, 0);
			
				__m128d poti=_mm_setr_pd(0, 0);
				
				i==i1_ ? _j1=j1_ : _j1=i+1;
				if ((_j1==i+1) && (i%2==0)) _j1=i;
				if (_j1<jL) _j1=jL;
		
				i==i2_ ? _j2=j2_ : _j2=N-1;
				if (jR<_j2) _j2=jR;
				
				if (_j1>_j2) continue;
				if (_j2%2==0) _j2++;
					
				// to spread j tasks among nn2 and nn part
				if (_j1<n-1) {
					
					_j11=_j1; 
					if (n-1<_j2){ _j12=n-1; _j21=n; _j22=_j2; } 
							 else	{ _j21=_j22=0; _j12=_j2;      }	
				
				}  else { 
					
					_j11=_j12=0;	
					_j21=_j1; _j22=_j2; 
					
				}
				
				// nn2
				for (uint j=_j11; j<_j12; j+=2) {
							
					__m128d x0j=_mm_load_pd(&x0[j]); 
					__m128d dx0ij=_mm_sub_pd(x0i, x0j);
					__m128d sqrdxij0=_mm_mul_pd(dx0ij, dx0ij);
					
					__m128d x1j=_mm_load_pd(&x1[j]);
					__m128d dx1ij=_mm_sub_pd(x1i, x1j);
					__m128d sqrdxij1=_mm_mul_pd(dx1ij, dx1ij);
					
					__m128d x2j=_mm_load_pd(&x2[j]);	
					__m128d dx2ij=_mm_sub_pd(x2i, x2j);
 					__m128d sqrdxij2=_mm_mul_pd(dx2ij, dx2ij);

					
 					__m128d v0j=_mm_load_pd(&v0[j]); 
 					__m128d dv0ij=_mm_sub_pd(v0i, v0j);
 					
					__m128d v1j=_mm_load_pd(&v1[j]);
					__m128d dv1ij=_mm_sub_pd(v1i, v1j);
					
					__m128d v2j=_mm_load_pd(&v2[j]);	
					__m128d dv2ij=_mm_sub_pd(v2i, v2j);
					
					
					__m128d rij2 = _mm_add_pd(sqrdxij0, sqrdxij1);
					rij2 = _mm_add_pd(rij2, sqrdxij2);
					rij2 = _mm_add_pd(rij2, _eps2);
			
					__m128d rij = _mm_sqrt_pd(rij2);
					__m128d rij3 = _mm_mul_pd(rij2, rij);

					__m128d mj=_mm_load_pd(&m[j]);
					__m128d mimj = _mm_mul_pd(mi, mj);
					__m128d mimjrij3 = _mm_div_pd(mimj, rij3);
					
					__m128d dxijdvij0 =_mm_mul_pd(dx0ij, dv0ij);
					__m128d dxijdvij1 =_mm_mul_pd(dx1ij, dv1ij);
					__m128d dxijdvij2 =_mm_mul_pd(dx2ij, dv2ij);
					
					__m128d 
					xv3rij2 = _mm_add_pd(dxijdvij0, dxijdvij1);
					xv3rij2 = _mm_add_pd(xv3rij2, dxijdvij2);
					xv3rij2 = _mm_mul_pd(xv3rij2, _3);
					xv3rij2 = _mm_div_pd(xv3rij2, rij2);
					
					// POT 
					
					__m128d potij = _mm_mul_pd(mimjrij3, rij2);
 					poti =_mm_sub_pd(poti, potij);	
					__m128d potj =_mm_load_pd(&pot[j]);
					potj =_mm_sub_pd(potj, potij);
					_mm_store_pd(&pot[j], potj);
					
					
					__m128d f0ij = _mm_mul_pd(mimjrij3, dx0ij);
					f0i =_mm_sub_pd(f0i, f0ij);
					__m128d f0j=_mm_load_pd(&f0[j]);
					f0j =_mm_add_pd(f0j, f0ij);
					_mm_store_pd(&f0[j], f0j);
					
					__m128d f1ij = _mm_mul_pd(mimjrij3, dx1ij);
					f1i =_mm_sub_pd(f1i, f1ij);
					__m128d f1j=_mm_load_pd(&f1[j]);
					f1j =_mm_add_pd(f1j, f1ij);
					_mm_store_pd(&f1[j], f1j);
					
					__m128d f2ij = _mm_mul_pd(mimjrij3, dx2ij);
					f2i =_mm_sub_pd(f2i, f2ij);
					__m128d f2j=_mm_load_pd(&f2[j]);
					f2j =_mm_add_pd(f2j, f2ij);
					_mm_store_pd(&f2[j], f2j);



					
					__m128d tmp0 =	_mm_mul_pd(xv3rij2, dx0ij);
					tmp0 =	_mm_sub_pd(dv0ij, tmp0);
					__m128d j0ij = _mm_mul_pd(mimjrij3, tmp0);
					j0i =_mm_sub_pd(j0i, j0ij);
					__m128d j0j=_mm_load_pd(&j0[j]);
					j0j =_mm_add_pd(j0j, j0ij);
					_mm_store_pd(&j0[j], j0j);
					
					__m128d tmp1 =	_mm_mul_pd(xv3rij2, dx1ij);
					tmp1 = _mm_sub_pd(dv1ij, tmp1);
					__m128d j1ij = _mm_mul_pd(mimjrij3, tmp1);
					j1i =_mm_sub_pd(j1i, j1ij);
					__m128d j1j=_mm_load_pd(&j1[j]);
					j1j =_mm_add_pd(j1j, j1ij);
					_mm_store_pd(&j1[j], j1j);
					
					__m128d tmp2 =	_mm_mul_pd(xv3rij2, dx2ij);
					tmp2 = _mm_sub_pd(dv2ij, tmp2);
					__m128d j2ij = _mm_mul_pd(mimjrij3, tmp2);
					j2i =_mm_sub_pd(j2i, j2ij);
					__m128d j2j=_mm_load_pd(&j2[j]);
					j2j =_mm_add_pd(j2j, j2ij);
					_mm_store_pd(&j2[j], j2j);
					
 					
				}

				#ifdef DEBUG
				printf("j2=[%u, %u]\n", _j21, _j22);
				#endif
				
				//nn
				for (uint j=_j21; j<_j22; j+=2) {
					
					
					
					__m128d x0j=_mm_load_pd(&x0[j]); 
					__m128d dx0ij=_mm_sub_pd(x0i, x0j);
					__m128d sqrdxij0=_mm_mul_pd(dx0ij, dx0ij);
					__m128d v0j=_mm_load_pd(&v0[j]); 
 					__m128d dv0ij=_mm_sub_pd(v0i, v0j);
 										
					__m128d x1j=_mm_load_pd(&x1[j]);
					__m128d dx1ij=_mm_sub_pd(x1i, x1j);
					__m128d sqrdxij1=_mm_mul_pd(dx1ij, dx1ij);
					__m128d v1j=_mm_load_pd(&v1[j]);
					__m128d dv1ij=_mm_sub_pd(v1i, v1j);
					
					
					__m128d x2j=_mm_load_pd(&x2[j]);	
					__m128d dx2ij=_mm_sub_pd(x2i, x2j);
 					__m128d sqrdxij2=_mm_mul_pd(dx2ij, dx2ij);
					__m128d v2j=_mm_load_pd(&v2[j]);	
					__m128d dv2ij=_mm_sub_pd(v2i, v2j);
					
			
					__m128d rij2 = _mm_add_pd(sqrdxij0, sqrdxij1);
					rij2 = _mm_add_pd(rij2, sqrdxij2);
					rij2 = _mm_add_pd(rij2, _eps2);
			
					__m128d rij = _mm_sqrt_pd(rij2);
					__m128d rij3 = _mm_mul_pd(rij2, rij);

					__m128d mj=_mm_load_pd(&m[j]);
					__m128d mimj = _mm_mul_pd(mi, mj);
					__m128d mimjrij3 = _mm_div_pd(mimj, rij3);
					
					// POT 
					
					__m128d potij = _mm_mul_pd(mimjrij3, rij2);
					poti =_mm_sub_pd(poti, potij);	
					
					
					__m128d f0ij = _mm_mul_pd(mimjrij3, dx0ij);
					f0i =_mm_sub_pd(f0i, f0ij);
					__m128d dxijdvij0 =_mm_mul_pd(dx0ij, dv0ij);
					
					__m128d f1ij = _mm_mul_pd(mimjrij3, dx1ij);
					f1i =_mm_sub_pd(f1i, f1ij);
					__m128d dxijdvij1 =_mm_mul_pd(dx1ij, dv1ij);
					
					__m128d f2ij = _mm_mul_pd(mimjrij3, dx2ij);
					f2i =_mm_sub_pd(f2i, f2ij);
					__m128d dxijdvij2 =_mm_mul_pd(dx2ij, dv2ij);
					
					__m128d 
					xv3rij2 = _mm_add_pd(dxijdvij0, dxijdvij1);
					xv3rij2 = _mm_add_pd(xv3rij2, dxijdvij2);
					xv3rij2 = _mm_mul_pd(xv3rij2, _3);
					xv3rij2 = _mm_div_pd(xv3rij2, rij2);
					
					__m128d 
					tmp0 = _mm_mul_pd(xv3rij2, dx0ij);
					tmp0 = _mm_sub_pd(dv0ij, tmp0);
					__m128d j0ij = _mm_mul_pd(mimjrij3, tmp0);
					j0i =_mm_sub_pd(j0i, j0ij);
					
					__m128d 
					tmp1 = _mm_mul_pd(xv3rij2, dx1ij);
					tmp1 = _mm_sub_pd(dv1ij, tmp1);
					__m128d j1ij = _mm_mul_pd(mimjrij3, tmp1);
					j1i =_mm_sub_pd(j1i, j1ij);
					
					__m128d 
					tmp2 = _mm_mul_pd(xv3rij2, dx2ij);
					tmp2 = _mm_sub_pd(dv2ij, tmp2);
					__m128d j2ij = _mm_mul_pd(mimjrij3, tmp2);
					j2i =_mm_sub_pd(j2i, j2ij);
					
				}

				_mm_store_pd(tmp, f0i);
				f0[i]+=tmp[0]+tmp[1];
		
				_mm_store_pd(tmp, f1i);
				f1[i]+=tmp[0]+tmp[1];
		
				_mm_store_pd(tmp, f2i);
				f2[i]+=tmp[0]+tmp[1];
				
				
				_mm_store_pd(tmp, j0i);
				j0[i]+=tmp[0]+tmp[1];
		
				_mm_store_pd(tmp, j1i);
				j1[i]+=tmp[0]+tmp[1];
		
				_mm_store_pd(tmp, j2i);
				j2[i]+=tmp[0]+tmp[1];

				
				
				_mm_store_pd(tmp, poti);
				if (i==_j11) tmp[0]+=2*SQR(m[i])/sqrt(eps2);	// i==j problem :: tempory :: substitute when will separate potential....
				pot[i]+=tmp[0]+tmp[1];
				
			}
		}
		
	}	
};

}
#endif